Para el análisis se utilizarán los datos abiertos proporcionados por el portal Portal de datos abiertos del Ayuntamiento de Madrid. Concretamente, los facilitados por el Sistema Integral de la Calidad del Aire del Ayuntamiento de Madrid, que pone a disposición los datos de los contaminantes registrados por las estaciones de medición situadas en Madrid. Los datos son de frecuencia horaria por anualidades desde 2001 y los datos se actualizan de forma mensual.
Además, existe una API para obtener los datos en tiempo real.
NOTA IMPORTANTE:
los conjuntos de datos de la calidad de aire son complejos y en algunos casos los datos no pueden utilizarse tal cual y pueden requerir consideraciones cuidadosas antes de llegar a cualquier conclusión. Debe prestarse atención a la existencia de subgrupos.
¿Dónde se han medido los contaminantes del aire?
Aunque los datos en bruto no son inmediatos de tratar, existe una amplia comunidad de personas que comparten su código y tratan los datos de forma que sean más inmediatos para su uso.
Existe código en GitHub como el repositorio michal0091/aire_madrid que trata los datos brutos, facilitando su uso. Las ventajas del uso del código compartido de forma pública son:
la transparencia: tenemos acceso al código y sus versiones
mantenimiento: en GitHUb es muy fácil saber cuándo fue la última actualización
interacción: si detectamos un error siempre se puede abrir un issue para avisar al propietario o incluso solicitar un permiso de push y corregir el error
Cargamos o instalamos las librerías si no las tenemos instaldas
# Establecemos el directorio de trabajo en la misma carpeta del repo
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
# Uncoment to install
# install.packages("raster", dependencies = TRUE) # no permitir el binario
# install.packages("terra", dependencies = TRUE)
# install.packages("gstat", dependencies = TRUE)
#
# install.packages(c("tidyverse", "lubridate", "ggridges", "ggrepel", "gganimate", "viridis", "hrbrthemes", "ggstatsplot", "sf", "mapSpain"), dependencies = TRUE)
# Libraries
library(gstat) # Estadística espacial
library(tidyverse) # Kit de libraries para data science
library(lubridate) # Manejo de fechas
# Libaries para gráficos
library(ggridges)
library(ggrepel)
library(gganimate)
library(viridis)
library(hrbrthemes)
library(ggstatsplot)
# Libraries para manejo de datos espaciales
library(sf)
library(mapSpain)
library(raster)
# Data
air_mad <- readr::read_rds("dt_daily_mean_2011.RDS")
Una buena manera para ver los datos es presentar sus primeras líneas:
head(air_mad)
## estaciones id id_name longitud latitud nom_mag nom_abv
## 1 E08: Escuelas Aguirre E08 Escuelas Aguirre -3.682316 40.42155 Benceno BEN
## 2 E08: Escuelas Aguirre E08 Escuelas Aguirre -3.682316 40.42155 Benceno BEN
## 3 E08: Escuelas Aguirre E08 Escuelas Aguirre -3.682316 40.42155 Benceno BEN
## 4 E08: Escuelas Aguirre E08 Escuelas Aguirre -3.682316 40.42155 Benceno BEN
## 5 E08: Escuelas Aguirre E08 Escuelas Aguirre -3.682316 40.42155 Benceno BEN
## 6 E08: Escuelas Aguirre E08 Escuelas Aguirre -3.682316 40.42155 Benceno BEN
## ud_med fecha daily_mean zona tipo
## 1 µg/m3 2011-01-01 0.9375000 Interior M30 Tráfico
## 2 µg/m3 2011-01-02 1.0666667 Interior M30 Tráfico
## 3 µg/m3 2011-01-03 1.6666667 Interior M30 Tráfico
## 4 µg/m3 2011-01-04 1.1416667 Interior M30 Tráfico
## 5 µg/m3 2011-01-05 1.0416667 Interior M30 Tráfico
## 6 µg/m3 2011-01-06 0.4166667 Interior M30 Tráfico
Para ver la dimensión de la tabla de datos:
dim(air_mad)
## [1] 521388 12
Luego se trata de 521388 filas y 12 columnas.
Una de las funciones más básicas para ver la estructura de los datos
es la función str(). Éste es el que debería de ser uno de
los primeros pasos al analizar un dataset nuevo:
str(air_mad)
## Classes 'data.table' and 'data.frame': 521388 obs. of 12 variables:
## $ estaciones: chr "E08: Escuelas Aguirre" "E08: Escuelas Aguirre" "E08: Escuelas Aguirre" "E08: Escuelas Aguirre" ...
## $ id : chr "E08" "E08" "E08" "E08" ...
## $ id_name : chr "Escuelas Aguirre" "Escuelas Aguirre" "Escuelas Aguirre" "Escuelas Aguirre" ...
## $ longitud : num -3.68 -3.68 -3.68 -3.68 -3.68 ...
## $ latitud : num 40.4 40.4 40.4 40.4 40.4 ...
## $ nom_mag : chr "Benceno" "Benceno" "Benceno" "Benceno" ...
## $ nom_abv : chr "BEN" "BEN" "BEN" "BEN" ...
## $ ud_med : chr "µg/m3" "µg/m3" "µg/m3" "µg/m3" ...
## $ fecha : Date, format: "2011-01-01" "2011-01-02" ...
## $ daily_mean: num 0.938 1.067 1.667 1.142 1.042 ...
## $ zona : chr "Interior M30" "Interior M30" "Interior M30" "Interior M30" ...
## $ tipo : chr "Tráfico" "Tráfico" "Tráfico" "Tráfico" ...
## - attr(*, "sorted")= chr [1:9] "estaciones" "id" "id_name" "longitud" ...
## - attr(*, ".internal.selfref")=<externalptr>
En la salida observamos que se trata de un data.frame del tipo data.table, la dimensión de los datos, las variables y su clase. Las clases que aparecen son:
chr: "character" se trata de una
cadena de texto
num: "numeric" se trata de un
vector numérico
Date: "Date" formato de fecha.
Internamente, los objetos Date se almacenan como el número
de días desde el 1 de enero de 1970, utilizando números negativos para
fechas anteriores. La función as.numeric() puede utilizarse
para convertir un objeto Date a su forma interna.
Dado que queremos trabajar con tidyverse vamos a convertir este tipo de tabla en una tibble.
air_mad <- as_tibble(air_mad)
Veamos si la clase ha cambiado:
class(air_mad)
## [1] "tbl_df" "tbl" "data.frame"
Ahora estamos trabajando con un objeto tibble.
La función summary() proporciona una salida estadística
estándar para cada columna:
summary(air_mad)
## estaciones id id_name longitud
## Length:521388 Length:521388 Length:521388 Min. :-3.775
## Class :character Class :character Class :character 1st Qu.:-3.713
## Mode :character Mode :character Mode :character Median :-3.689
## Mean :-3.683
## 3rd Qu.:-3.661
## Max. :-3.580
##
## latitud nom_mag nom_abv ud_med
## Min. :40.35 Length:521388 Length:521388 Length:521388
## 1st Qu.:40.40 Class :character Class :character Class :character
## Median :40.42 Mode :character Mode :character Mode :character
## Mean :40.43
## 3rd Qu.:40.46
## Max. :40.52
##
## fecha daily_mean zona tipo
## Min. :2011-01-01 Min. : 0.000 Length:521388 Length:521388
## 1st Qu.:2013-10-31 1st Qu.: 2.792 Class :character Class :character
## Median :2016-08-30 Median : 13.208 Mode :character Mode :character
## Mean :2016-08-30 Mean : 26.065
## 3rd Qu.:2019-07-01 3rd Qu.: 34.125
## Max. :2022-04-30 Max. :4063.458
## NA's :15615
Otra forma muy común en ciencia de datos de hacer un summary es
mediante la función skim() de la librería
skmimr:
skimr::skim(air_mad)
| Name | air_mad |
| Number of rows | 521388 |
| Number of columns | 12 |
| _______________________ | |
| Column type frequency: | |
| character | 8 |
| Date | 1 |
| numeric | 3 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| estaciones | 0 | 1 | 11 | 26 | 0 | 23 | 0 |
| id | 0 | 1 | 3 | 3 | 0 | 23 | 0 |
| id_name | 0 | 1 | 6 | 21 | 0 | 23 | 0 |
| nom_mag | 0 | 1 | 7 | 21 | 0 | 10 | 0 |
| nom_abv | 0 | 1 | 2 | 5 | 0 | 10 | 0 |
| ud_med | 0 | 1 | 5 | 5 | 0 | 1 | 0 |
| zona | 0 | 1 | 7 | 12 | 0 | 5 | 0 |
| tipo | 0 | 1 | 5 | 10 | 0 | 3 | 0 |
Variable type: Date
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| fecha | 0 | 1 | 2011-01-01 | 2022-04-30 | 2016-08-30 | 4138 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| longitud | 0 | 1.00 | -3.68 | 0.05 | -3.77 | -3.71 | -3.69 | -3.66 | -3.58 | ▂▆▇▂▃ |
| latitud | 0 | 1.00 | 40.43 | 0.04 | 40.35 | 40.40 | 40.42 | 40.46 | 40.52 | ▁▆▇▅▂ |
| daily_mean | 15615 | 0.97 | 26.07 | 40.23 | 0.00 | 2.79 | 13.21 | 34.12 | 4063.46 | ▇▁▁▁▁ |
NOTA IMPORTANTE:
Además, existen paquetes como DataExplorer y
dlookr que generan informes automáticos con los principales
descriptivos.
Con las funciones del Tidyverse represente la evolución de todos los contaminantes medidos por las estaciones de monitoreo de la ciudad de Madrid en el periodo (2011-2022).
plot_air_mad <- air_mad %>%
group_by(fecha, nom_mag) %>%
summarise(media_estaciones = mean(daily_mean, na.rm = TRUE)) %>%
ggplot(aes(x = fecha, y = media_estaciones)) +
geom_line() +
geom_smooth() +
facet_wrap( ~ nom_mag, scales = "free_y", ncol = 2)
plot_air_mad
NOTA IMPORTANTE:
La representación gráfica es una de las herramientas más poderosas en el EDA siempre y cuando esté bien definida.
Con las funciones del Tidyverse y la librería lubridate
represente la evolución de todos los contaminantes medidos por las
estaciones de monitoreo de la ciudad de Madrid en el periodo (2011-2022)
y personalice los plots con distintos colores.
air_mad %>%
group_by(semana=floor_date(fecha,unit = "week"), nom_mag) %>%
summarise(media_estaciones=mean(daily_mean, na.rm=TRUE)) %>%
ggplot(aes(x=semana, y=media_estaciones))+
geom_line(aes(color=nom_mag))+
geom_smooth(size=0.5, color="black")+
scale_color_brewer(palette = "Paired")+
labs(x=NULL, y="(µg/m3)", title = "Evolución de partículas contaminantes en Madrid",
subtitle = "Concentración media semanal en las estaciones de medición",
caption = "Fuente: Portal de datos abiertos del Ayuntamiento de Madrid" )+
theme_minimal()+
theme(legend.position = "none")+
facet_wrap(~nom_mag, scales = "free_y", ncol = 2)
Antes de continuar. No nos olvidemos de los datos faltantes, a veces un importante proble en ciencia de datos… ¿Existe algún patrón en los NAs?
¿Cuántos NAs tengo en mis datos?
sum(is.na(air_mad))
## [1] 15615
¿Dónde están los NAs en mis datos?
na_table <- air_mad %>%
mutate(isna = is.na(daily_mean)) %>%
ggplot(aes(x = fecha, y = id_name, fill = isna)) +
geom_raster(alpha = 0.8) +
theme(legend.position = "bottom")
na_table
¿Cuándo se han producido los Nas?
na_fechas <- air_mad %>%
filter(is.na(daily_mean)) %>%
group_by(id_name, fecha = floor_date(fecha, unit = "month")) %>%
summarise(num_nas = sum(is.na(daily_mean))) %>%
ungroup() %>%
mutate(fecha = paste0(year(fecha), "-", month(fecha, label = TRUE, abbr = FALSE))) %>%
arrange(desc(num_nas))
Opción 1: Eliminar los NAs
air_mad_clean <- air_mad %>%
drop_na()
summary(is.na(air_mad_clean))
## estaciones id id_name longitud
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:505773 FALSE:505773 FALSE:505773 FALSE:505773
## latitud nom_mag nom_abv ud_med
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:505773 FALSE:505773 FALSE:505773 FALSE:505773
## fecha daily_mean zona tipo
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:505773 FALSE:505773 FALSE:505773 FALSE:505773
Opción 2: Imputar NAs con día anterior
air_mad_dia_antes <- air_mad %>%
arrange(estaciones, nom_abv, fecha) %>%
fill(daily_mean)
summary(is.na(air_mad_dia_antes))
## estaciones id id_name longitud
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:521388 FALSE:521388 FALSE:521388 FALSE:521388
## latitud nom_mag nom_abv ud_med
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:521388 FALSE:521388 FALSE:521388 FALSE:521388
## fecha daily_mean zona tipo
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:521388 FALSE:521388 FALSE:521388 FALSE:521388
A lo largo de los años, como se ha podido apreciar, el nivel de todos los contaminantes estudiados ha descendido. Sin embargo, hay dos que aun presentan elevado número de superaciones del estandar legal y que son muy dañinos para la saludo: PM10 y NOx.
Una herramienta muy útil para tener una visión general de estos contaminantes viendo el calendario como un heatmap: calendar heatmap
calendar_plot <- air_mad %>%
group_by(fecha, nom_mag, nom_abv, ud_med) %>%
summarize(valor_promedio = mean(daily_mean, na.rm = T))
# Dates as factors
months <-
seq.Date(from = as.Date("2022-01-01"),
length.out = 12,
by = "month") %>% format("%B")
wdays <-
seq.Date(from = as.Date("2022-05-30"),
length.out = 7,
by = "day") %>% format("%A")
calendar_plot <- calendar_plot %>%
mutate(year = format(fecha, "%Y"),
month = factor(format(fecha, "%B"), levels = months, labels = months),
wday = factor(weekdays(fecha), levels = wdays, labels = wdays),
week = as.numeric(format(fecha, "%W")))
calendar_plot <- calendar_plot %>%
group_by(year, month) %>%
mutate(wmonth = 1 + week - min(week))
Calendar heatmap para NOx
i_mag <- "Óxidos de Nitrógeno"
fill_title <- calendar_plot %>%
filter(nom_mag == i_mag & year >= 2011) %>%
ungroup() %>%
distinct(paste(unique(nom_abv), unique(ud_med)))
calendar_plot %>%
filter(nom_mag == i_mag & year >= 2011) %>%
ggplot(aes(
x = wmonth,
y = reorder(wday, -as.numeric(wday)),
fill = valor_promedio
)) +
geom_tile(colour = "white") +
facet_grid(year ~ month) +
scale_fill_gradient(low = "yellow", high = "red", ) +
scale_x_continuous(breaks = 1:5, limits = c(0, 6)) +
labs(
x = "Semana del mes",
y = NULL,
title = paste0("Concentración de ", i_mag, " por día de la semana"),
fill = fill_title,
caption = "Fuente: Red de Vigilancia de la Calidad del Aire del Ayto. de Madrid"
)
Calendar heatmap para PM10
i_mag <- "Partículas < 10 µm"
fill_title <- calendar_plot %>%
filter(nom_mag == i_mag & year >= 2011) %>%
ungroup() %>%
distinct(paste(unique(nom_abv), unique(ud_med)))
calendar_plot %>%
filter(nom_mag == i_mag & year >= 2011) %>%
ggplot(aes(
x = wmonth,
y = reorder(wday, -as.numeric(wday)),
fill = valor_promedio
)) +
geom_tile(colour = "white") +
facet_grid(year ~ month) +
scale_fill_gradient(low = "yellow", high = "red", ) +
scale_x_continuous(breaks = 1:5, limits = c(0, 6)) +
labs(
x = "Semana del mes",
y = NULL,
title = paste0("Concentración de ", i_mag, " por día de la semana"),
fill = fill_title,
caption = "Fuente: Red de Vigilancia de la Calidad del Aire del Ayto. de Madrid"
)
Seleccione los dos contaminantes más problemáticos en la ciudad de Madrid (PM10 y NOx) y representelos en forma de serie temporal.
air_mad_pm10_nox <- air_mad %>%
filter(nom_abv %in% c("PM10", "NOx"))
plot_pm10_nox <- air_mad_pm10_nox %>%
group_by(fecha, nom_mag) %>%
summarise(media_estaciones = mean(daily_mean, na.rm = TRUE)) %>%
ggplot(aes(x = fecha, y = media_estaciones)) +
geom_line(aes(color = nom_mag)) +
geom_smooth(size = 0.5, color = "black", se = FALSE) +
scale_color_brewer(palette = "Paired") +
labs(
x = NULL,
y = "(µg/m3)",
title = "Evolución semanal de partículas contaminantes (PM10 y NOx) en Madrid",
subtitle = "Concentración media semanal en las estaciones de medición",
caption = "Fuente: Portal de datos abiertos del Ayuntamiento de Madrid"
) +
theme_minimal() +
theme(legend.position = "none") +
facet_wrap( ~ nom_mag, scales = "free_y", ncol = 1)
plot_pm10_nox
NOTA IMPORTANTE:
La función ggplotly() de la librería plotly
permite hacer fácilmente gráficos interactivos.
¿Por qué no hacer el anterior gráfico interactivo con
plotly?
plotly:: ggplotly(plot_pm10_nox)
Una vez tratados los aspectos temporales de la variable, ¿por qué no analizar la dimensión espacial y el tipo de estación de monitoreo?
Gráfico de densidad (Ridgeline) PM10 por tipo de
estación
air_mad %>%
filter(nom_abv == "PM10") %>%
ggplot(aes(x = daily_mean, y = tipo, fill = tipo)) +
geom_density_ridges() +
theme_ridges() +
scale_x_continuous(limits = c(0, 100), breaks = seq(0, 100, 10)) +
theme(legend.position = "none")
Gráfico de densidad (Ridgeline) PM10 por zona de
calidad del aire
air_mad %>%
filter(nom_abv == "NOx") %>%
ggplot(aes(x = daily_mean, y = zona, fill = ..x..)) +
geom_density_ridges_gradient(scale = 2, rel_min_height = 0.001) +
scale_fill_viridis(option = "C") +
scale_x_continuous(limits = c(0, 75), ) +
labs(title = 'Concentración de PM10 µm por tipo de estación en Madrid (2011-2022)') +
xlab("Concentración diaria (µg/m3)") +
theme_ipsum() +
theme(legend.position = "none",
strip.text.x = element_text(size = 8))
Repita el análisis anterior para el PM10 y compruebe si la zona “Interior M-30” es la que presenta mayor contaminación o no
¿Y por zona de calidad del aire?
air_mad %>%
filter(nom_abv == "NOx") %>%
ggplot(aes(x = daily_mean, y = zona, fill = ..x..)) +
geom_density_ridges_gradient(scale = 2, rel_min_height = 0.001) +
scale_fill_viridis(option = "C") +
scale_x_continuous(limits = c(0, 75), ) +
labs(title = 'Concentración de NOx µm por zona en Madrid (2011-2022)') +
xlab("Concentración diaria (µg/m3)") +
theme_ipsum() +
theme(legend.position = "none",
strip.text.x = element_text(size = 8))
Compruebe con un Análisis de la Varianza si los niveles de concentración de PM10 dependen de las variables tipo y zona
air_mad_anova <- air_mad %>%
filter(nom_abv == "NOx") %>%
drop_na()
anova <- aov(data = air_mad_anova, daily_mean ~ zona + tipo)
summary(anova)
## Df Sum Sq Mean Sq F value Pr(>F)
## zona 4 22508217 5627054 1540.1 <2e-16 ***
## tipo 2 5489043 2744521 751.2 <2e-16 ***
## Residuals 94975 347009104 3654
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova2 <- lm(data = air_mad_anova, daily_mean ~ zona + tipo)
summary(anova2)
##
## Call:
## lm(formula = daily_mean ~ zona + tipo, data = air_mad_anova)
##
## Residuals:
## Min 1Q Median 3Q Max
## -91.32 -36.10 -16.78 13.51 679.90
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 61.4511 0.4653 132.070 < 2e-16 ***
## zonaNoreste -1.8855 0.6270 -3.007 0.00264 **
## zonaNoroeste -11.5701 1.3126 -8.814 < 2e-16 ***
## zonaSureste -2.5472 0.6504 -3.916 9e-05 ***
## zonaSuroeste 23.2637 0.6504 35.767 < 2e-16 ***
## tipoSuburbanas -20.2272 1.0315 -19.609 < 2e-16 ***
## tipoTráfico 17.2328 0.5154 33.434 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 60.45 on 94975 degrees of freedom
## Multiple R-squared: 0.07466, Adjusted R-squared: 0.0746
## F-statistic: 1277 on 6 and 94975 DF, p-value: < 2.2e-16
# select(PM10) by_group (estación de monitoreo) y gráfico de violín (ggstatsplot) para ver en el periodo cual ha sido la estación que ha registrado valores más altos.
air_mad %>%
filter(nom_abv == "PM10") %>%
filter(daily_mean < 1000) %>%
filter(between(
fecha,
left = as.Date("2021-05-01"),
right = as.Date("2022-04-30")
)) %>%
ggbetweenstats(
x = id_name,
y = daily_mean,
xlab = "Estación de medida",
ylab = "Concentración diaria (µg/m3)",
plot.type = "boxviolin",
# outlier.tagging = TRUE,
# outlier.coef = 1.5,
# outlier.label = fecha,
# outlier.label.args = list(color = "red", size=1),
title = "Grafico comparativo entre estaciones de concentración de PM10",
subtitle = "Madrid. Mayo 2021 - abril 2022.",
caption = "Fuente: Portal de datos abiertos del Ayuntamiento de Madrid",
)
¿Qué pasó en la semana de la calima con el PM10?
particulas <- c("Partículas < 2.5 µm", "Partículas < 10 µm")
calima <- air_mad %>%
filter(nom_mag %in% particulas &
fecha %in% seq.Date(as.Date("2022-03-01"), by = "day", length.out = 31)) %>%
group_by(fecha, id, id_name, nom_mag, nom_abv, ud_med) %>%
summarize(valor_promedio = mean(daily_mean, na.rm = T))
max_2.5 <- calima %>%
ungroup() %>%
filter(nom_mag == particulas[1]) %>%
slice(which.max(valor_promedio))
max_10 <- calima %>%
ungroup() %>%
filter(nom_mag == particulas[2]) %>%
slice(which.max(valor_promedio))
calima %>%
ggplot(aes(fecha, valor_promedio, colour = nom_mag)) +
geom_jitter() +
geom_smooth(
method = "loess",
span = .5,
se = FALSE,
show.legend = FALSE
) +
scale_x_date(breaks = seq.Date(as.Date("2022-03-01"), by = "week", length.out = 5),
date_labels = "%d-%b") +
scale_color_manual(values = c("#261606", "#DD9C4A")) +
geom_label_repel(
data = max_10,
mapping = aes(fecha, valor_promedio, label = paste(id_name)),
show.legend = FALSE
) +
geom_label_repel(
data = max_2.5,
mapping = aes(fecha, valor_promedio, label = paste(id_name)),
show.legend = FALSE
) +
labs(
title = "Registro de partículas durante el mes de marzo 2022",
subtitle = "Madrid",
x = NULL,
y = unique(calima$ud_med),
color = NULL,
caption = "Fuente: Red de Vigilancia de la Calidad del Aire del Ayto. de Madrid\nElaboración: Michal Kinel"
) +
theme(
legend.position = 'bottom',
panel.background = element_rect(
fill = "#f4dbb3",
colour = "#f4dbb3",
size = 0.5,
linetype = "solid"
),
plot.background = element_rect(fill = "#DD9C4A"),
text = element_text(family = "helvetica", colour = "#261606"),
legend.background = element_rect(fill = "#DD9C4A"),
legend.key = element_rect(fill = "#f4dbb3", color = NA)
)
# idw
mad_sf <- esp_get_munic(munic = "^Madrid$", epsg = 4326)
marzo_pm10 <- air_mad %>%
filter(nom_abv == "PM10" &
fecha >= as.Date("2022-03-13") & fecha <= as.Date("2022-03-17")) %>%
drop_na()
madrid_estaciones_sf <- st_as_sf(marzo_pm10,
coords = c("longitud", "latitud"),
crs = 4326)
ggplot(madrid_estaciones_sf) +
geom_sf(data = mad_sf,
fill = "grey95") +
geom_sf(
aes(fill = daily_mean),
shape = 21,
size = 5,
alpha = .7
) +
labs(fill = "PM10") +
scale_fill_viridis_c() +
theme_void() +
labs(title = "PM10: {current_frame}") +
transition_manual(fecha) +
ease_aes('linear') +
theme(
plot.title = element_text(size = 12,
face = "bold"),
plot.subtitle = element_text(size = 8,
face = "italic")
)
# We choose to project our objects to ETRS89 / UTM zone 30N EPSG:25830, which provides projected x and y values in meters and maximizes the accuracy for Spain.
madrid_estaciones_utm <- st_transform(madrid_estaciones_sf, 25830)
mad_utm <- st_transform(mad_sf, 25830)
test_day <- madrid_estaciones_utm %>% filter(fecha == "2022-03-13")
extent_mad_utm <- extent(mad_utm)
grid_mad_utm <-
expand.grid(
x = seq(
from = round(extent_mad_utm@xmin),
to = round(extent_mad_utm@xmax),
by = 500
),
y = seq(
from = round(extent_mad_utm@ymin),
to = round(extent_mad_utm@ymax),
by = 500
)
)
coordinates(grid_mad_utm) <- ~ x + y
grid_mad_utm <- st_as_sf(grid_mad_utm)
st_crs(grid_mad_utm) <- st_crs(mad_utm)
grid_mad_utm <- as(grid_mad_utm, "Spatial")
gridded(grid_mad_utm) <- TRUE
neighbors <- length(test_day)
beta <- 2
interp_pm10 = gstat::gstat(
formula = daily_mean ~ 1,
# intercept only model
data = test_day,
nmax = neighbors,
set = list(idp = beta)
)
grid_interp_pm10 <-
predict(object = interp_pm10, newdata = grid_mad_utm)
## [inverse distance weighted interpolation]
plot(grid_interp_pm10, main = "IDW PM10 µg/m³")
plot(st_geometry(mad_utm), add = TRUE, border = "white")
plot(
st_geometry(madrid_estaciones_utm),
add = TRUE,
pch = 19,
cex = 0.5,
col = "green"
)